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We show how standard Metadynamics coupled with classical Molecular Dynamics can be successfully ap¬ 
plied to sample the configurational and free energy space of metallic and bimetallic nanopclusters via the 
implementation of collective variables related to the pair distance distribution function of the nanoparticle 
itself. As paradigmatic examples we show an application of our methodology to Agi 47 , Pti 47 and their alloy 
Ags/ieZ/Ptcore ^t 1:1 and 2:1 chemical compositions. The proposed scheme is not only able to reproduce known 
structural transformation pathways, as the five and the six square-diamond mechanisms both in pure and 
core-shell nanoparticles but also to predict a new route connecting icosahedron to anti-cuboctahedron. 


I. INTRODUCTION 

Mono- and bi-metallic nanoparticles (mNPs) find a 
wide number of applications ranging from catalysis and 
biomedicine to optoelectronics and magnetic data stor¬ 
age due to their high surface to volume ratio, anisotropy 
and d-band shiftNanoclusters’ chemophysical prop¬ 
erties in fact strongly depend on the interplay between 
their size, morphology, and chemical composition. Un¬ 
derstanding the relative stability of a configuration de¬ 
pends not only on weather it is a local minimum on the 
potential energy landscape but also on complex entropic 
contribution difficult to address experimentall}!^. An es¬ 
timate of the magnitude of free energy barrier and in¬ 
sights on the mechanisms involved in solid-solid struc¬ 
tural rearrangements among different minima of the free 
energy landscape may shed light on how mNPs chemical 
features vary due to ageing and external factors such as 
temperature or pressure. 

A variety of theoretical attempts have been presented 
in the literature to extend classical thermodynamics con¬ 
cepts to fully address and include the contributions deriv¬ 
ing by the large presence of surface atoms.^® Two are the 
general approaches to obtain a quantitative sampling of 
the nanoclusters’ free energy surface (FES) via computa¬ 
tional techniques. Double ended searches are based upon 
the foreknowledge of the initial and final point of the 
transition and consist in analysing accurately the tran¬ 
sition networks between the two target configurations, 
notable examples are the string methocPI, the discrete 
path sampling^ and the nudged elastic bancP. Temper¬ 
ature accelerated and biased sampling techniques repre¬ 
sent instead an opposite approach: an initial configura¬ 
tion of the system is excited or perturbed and forced to 
visit new and initially unknown attraction basins. Adap¬ 
tive biasiM forc^^, umbrella samplin^^ and parallel 
temperingi^ are all renown techniques based upon this 
idea. Both methodologies have found a wide application 
in the analysis and study of structural transformations in 
nanoclusters.l^^Hi^ Double ended methods weakness lies 
in the curse of dimensionality which afflicts the potential 
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energy surface (PES) minima basins network such that 
these algorithms need a larger and larger amount of com¬ 
putational expenses with increasing mNP size. Instead 
perturbation methods have been commonly used to de¬ 
tect order-disorder transition while rarely solid-solid ones 
have been characterized or simulated by these. 

In this paper we will show for the first time how 
Metadynamics can be successfully employed to sample 
the configurational space of relatively large metallic and 
bimetallic nanoparticles at finite temperatures and how 
the free energy barrier for various minimum-to-minimum 
can be estimated. Metadynamics (MetaD) algorithm 
combines the ionic dynamics with a history dependent 
potential exploited to accelerating rare events and to re¬ 
construct the free energy surface projections into an order 
parameters, named collective variables (CVs), space.l^^ 
Eirst introduced in the biophysics community to eluci¬ 
date ligand binding and protein folding phenomengP^ it 
has been demonstrated in the literature that this method 
can be also applied for finite inorganic nanos ystem s such 
as semiconductor/qua ntum dots nanopart icle d^ ^ * ^^ , alkali 
halides nanostr ncture P^^ Lennar d-Jone^^^^ or metal¬ 
lic nanopart icleJ^^^ of a very small size. However, it 
has been never applied extensively and systematically 
for metallic nanoclusters larger than few tens of atoms 
nor to bimetallic cases. We device molecular dynam¬ 
ics (MD), based on the second-moment approximation of 
the tight-binding potential (TBSMA), coupled to MetaD 
as a cost effective and accurate technique. Our MetaD 
scheme, with CVs based on the pair distribution function, 
has shown to be able to reproduced the well known five 
and six Diamond-Square-Diamond (DSD)I^ mechanism 
and for the first time the interconversion of an icosahe¬ 
dron into a anti-cuboctahedron is reproduced validating 
Mackay’s theoretical speculation.^^ 


II. MODELS AND METHOD 

We consider monometallic (Ag,Pt) and bimetallic 
(AgPt) nanoparticles at 1:1 and 2:1 chemical compo¬ 
sitions at a selected size of 147 atoms. At this size, 
noble and quasi-noble metals adopt geometrical closed 
shell polyhedra such as icosahedron (Ih), decahedron 











2 


(Dh), truncated octahedron or cubo-octahedron (Co) 
and rarely hexagonal close-packed geometries (hep), as 
reported in Figure As shown in details by Paz-Borbon 
et al!^, AgPt has a small size mismatch, with a ten¬ 
dency to segregation and a core-shell chemical ordering 
with silver in the uppermost layer is usually preferred. 



FIG. 1. Closed shell polyhedra at 147 atoms. Top row, from 
left to right: icosahedron, decahedron, cuboctahedron. Bot¬ 
tom row: hexagonal closed packed and anticuboctahedron. 


We perform classical molecular dynamics simulations, 
coupled with MetaD, where a velocity-Verlet algorithm is 
used for solving the Newton’s equations with a time step 
of 5 fs. An Anderesen thermostat is applied. Atomic 
interactions are modelled within the second moment ap¬ 
proximation of the tight binding theor}!^ (TBSMA), 
where the potential Vtbsma{x) is the sum of atomic con¬ 
tribution, 

^TBSMA ^ (1) 


Here the sum is extended up to the number of atoms 
Uy within an appropriate cut-off distance; a and b refer 
to the chemical species of the two atoms, is the bulk 
nearest neighbour distance, Aab and ^ab are related to 
the cohesive energy, while Pab and Qab and their product 
determine the range of the repulsive and of the attractive 
part of the potential. All the four parameters have been 
fitted using bulk properties. Parametrization of Ag and 
Pt interaction follows what suggested in Ref.^^ while for 
AgPt parameters are fitted according to Rem 

MetaD enhanced sampling algorithm is based upon 
the construction of a history-developing biasing poten¬ 
tial, AH, evolving as the sum of gaussians of height uj 
and width a added every to time interval, 

^(^) = Vtbsma{x) + AH {S{x),t) (2) 

[g(x)-s(tq]^ 

AV {S{x)A) = ^ bje 2ct2 

t'=tG,2tG.. 





S{x) represents a set of collective variables and defines 
the order parameter space on which the added gaussians 
lie; s{t') is the instanteneous value of the collective vari¬ 
able. When the collective variables have sampled all the 
conformation’s space they become fully diffusive and free 
energy along this order parameters can be reconstructed 
as the negative of the meta-potential AV{S{x)A)- The 
efficiency and physical faithfulness of MetaD is strongly 
related to the choice of a sensible set of collective vari¬ 
ables. The chosen metric should be able to distinguish 
between the various configurations the system visited 
during a run and the CVs should also be representative of 
all slow degree of motions involved in a structural trans¬ 
formation. Furthermore their number should be limited 
in such a way that we have a low dimensionality space 
to explore avoiding undesired computational expenses. If 
one of the above conditions is not respected the system 
can be forced to visit high energy regions and/or attrac¬ 
tion basin of interest may be hidden.^ 

The research of a good order parameter able to dif¬ 
ferentiate among the comple x con figuration space of a 
cluster is, however, not triviaP^^. On the other hand, 
an almost complete information on NP morphology is en¬ 
coded in its pair distance distribution function (PDDF), 
depicted in Figure where the four isomers shown in 
Figure are considered. 

Hence we broaden the restriction on our CV set ac¬ 
cording to the following criterion: it must able to clearly 
distinguish out the above mentioned geometries as these 
are the main structures of interest in our research. Thus, 
with a low computational cost and only a small loss in 
accuracy due to a partial correlation of the two CVs, 
we introduce specifically tailored order parameters cor¬ 
responding to window function, IFF, on the PDDF con¬ 
structed via a sigmoid function: 



where r^j,ro,do are respectively the distance between 
atom i and atom j, vq the window width, always set 
to 0.05 of the bulk lattice constant, and do is the charac¬ 
teristic distance, do is set to be 1.354 and 3.4 of the bulk 
FCC lattice parameter respectively for the stacking fault 
number (SFN) and the maximum pair distance difference 
(MPDD), as shown in Figure]^ 

A qualitative argument can be provided to support our 
collective variable choice. The MPDD is chosen to be set 
where the difference between configurations in the pair 
distribution function has a maximum. In such a way 
we are confident of being able to discriminate at least 
the main structural topologies. The stacking fault num¬ 
ber instead is related to a characteristic hep peak in the 
radial distribution function. Physically it is a topologi¬ 
cal defect obtained due to the intersection of two planes 
with different symmetry orientation. At the nanoscale. 
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FIG. 2. Pair distance distribution function for icosahedral 
(green), decahedral (blue), cuboctahedral (red), and hep 
(black) motifs of a 147 atoms cluster. The selected windows 
around the stacking fault number (SFN), and the maximum 
pair distribution difference (MPDD), are highlighted in grey 
and orange shadowing, respectively. 


as in the bulk, phase transformations happen via stacking 
faults, as demonstared for the diamond-square-diamond 
mechanismsConsequently a bias on this CV should 
not constrain the free energy surface exploration to un¬ 
physical structural transition or highly energetic states. 

We resort to the use of these two window functions 
on specific lattice distances because they adopt different 
values for each of the geometries of interest in our re¬ 
search and are computationally very cheap, being two 
body terms their calculations scales as N‘^. This will 
enable us to easily extend our methodology to larger sys¬ 
tems. On the other hand, this set of CVs can be applied 
with a partial confidence to bimetallic nanoalloys. The 
main problem lies in the fact that the presence of dif¬ 
ferent chemical species is not treated explicitly by the 
earlier on presented window functions, hence structural 
transformations involving chemical reordering are rarely 
reproduced biasing on those coordinates. However, due 
to the fact that Ag and Pt have a small size mismatch, 
and because we focus our attention on the transformation 
of a core-shell ordering, SFN and MPDD still capture the 
geometrical structural transition pathways of interest. 


A. Sampling and reconstructing the FES of a metallic 
nanoparticle 

During the course of the simulations the cluster topol¬ 
ogy is monitored on-the-fly during a MetaD-MD simula¬ 
tion via common neighbour analysis (CNA)P^. The CNA 
determines the local environment of all the nearest neigh¬ 
bours atom pair and produces a number of signatures 
percentages that characterises the whole cluster. A CNA 
signature is made by three integers (r, s, t): r is the num¬ 
ber of common nearest neighbour, s the number of bonds 


between the r-common nearest neighbours and t is the 
longest chain among the s-bonds. The percentage of each 
signature corresponds to features of the considered NP. 
Different CNA signatures are able to distinguish whether 
a pair of atoms is in a bulk environment or on the sur¬ 
face, whether the bulk is crystallographic or not, and if 
the pair belongs to a symmetry axis. A fast NP geome¬ 
try taxonomy may follow looking at the (4,2,2) and the 
(5,5,5) signatures. Their local neighbourhood is depicted 
in Figure]^ The (4,2,2) varies between 39% (Ih) to 25% 
(Dh and hep) and zero (Co); while the (5,5,5) ranges be¬ 
tween 5.2 % (Ih), 0.9 % (Dh) and zero (Co and hep). 



FIG. 3. Gommon neighbourhood, in yellow, for the red- 
coloured atomic pair. Left panel shows a (422) signature, 
associated with twin boundaries, while the right panel de¬ 
picts a (555) signature, which is related to pairs lying of a 
5-fold symmetry axis. 


A significant change in the collective variables and 
in the CNA signatures is therefore distinctive of struc¬ 
tural transition, as in the paradigmatic example of an 
Ag 92 @Pt 55 run sketched in Figure where three differ¬ 
ent basins, Ih (orange), Co (blue) and Dh (pink), are 
explored. Furthermore structures close to any presumed 
phase change are quenched in order to identify precisely 
the initial and final shape of the transition on the poten¬ 
tial energy surface. The PES profile obtained after that 
fast quenching is reported in the lower panel of Figure 

As it may happen that the diffusive regime is not 
reached, the free energy reconstruction is not merely the 
inverse of the added MetaD potential, AV in Eq. To 
obtain an overestimate of the free energy barrier sepa¬ 
rating two minima, we use the above on-the-fly analysis 
in order to monitor and collect after how many Gaus¬ 
sian depositions a transition takes place. Thus, an upper 
bound estimate of the free energy profile separating the 
two minima can be calculated as the negative of the re¬ 
pulsive MetaD potential deposited on the initial basin 
and needed to visit the following conformational mini¬ 
mum. 

In order to validate our approach and to compare our 
sampled free energy landscape with the corresponding 
potential energy one, the activation energy on the PES 
is calculated by means of a DETP^ where the relaxed 
structures obtained from the fast quenching of the con¬ 
figurations visited during a MetaD run serve as the ini¬ 
tial and endpoint of the transformation pathway, respec¬ 
tively. 
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FIG. 4. AgsheiiPtcore MetaD-MD simulation outcome and 
analysis. From top to bottom: Total energy difference 
{AEtot — Etot — in eV) and CV values are shown 

in the first three panels. In the lowest three panels, on-the- 
fly CNA analysis combined with a fast quenching for repro¬ 
ducing the map of all the basin encountered (AEquench = 
Equench “ Reference energy corresponds to the one 

of a quenched Ih . A colour scheme is applied for re¬ 

marking all the independent basin/shape assumed: Co-basin 
in blue, Ih-basin in orange, Dh-basin in pink. 


III. RESULTS 


Gaussians 0.25 eV high are deposited with a period 
of 20ps. The Gaussian width a is 15 along the SFN 
dimension and 10 for the MPDD. Such values are based 
on the standard deviations of the CVs during unbiased 
molecular dynamics simulations at 300 K and are chosen 
to assure the nearly spherical shape of the energy wells 
and thus an optimal conformational flooding. All the 
following results are obtained at room temperature of 
300 K. 



FIG. 5. Initial, saddle and final configurations of DSD mech¬ 
anism for Ih into Co (top row) and Ih into Dh (bottom row) 
transformation. Multicolored atoms delimit a facet in the 
original Ih. Rotation axis is shown. 


The detected structural transformations in the Ag clus¬ 
ters are Dh into Ih, Co into Ih and backwards, Ih into aCo 
and backwards. In the case of the Pt cluster only trans¬ 
formations into Ih geometry from Dh and Co have been 
observed before the exploration of highly defected struc¬ 
tures. The general solid-solid structural transformation 
mechanism between the above mentioned geometries is 
the so-called diamond - square - diamond (DSD) process 
as shown in Figures and where the motion of each 
atom is clearly identified via the different coloring. It 
consists of a stretching and rotation of triangular facets 
into a square by a collective dislocation of atoms along 
the surfaces involved. This dislocation corresponds to a 
rotation of different angles according to the initial and fi¬ 
nal configuration. Conversely a squared facet can trans¬ 
form into two triangular facets by the opposite move¬ 
ment. Five parallelograms are involved in this collective 
rearrangement in the case of the Dh-Ih transformation, 
and six in the case of the Co-Ih transformation. Figure 

n 



FIG. 6. Ih transforming into aGo via DSD mechanism. Sad¬ 
dle point is pictured in between the two. Multicolored atoms 
delimit a facet in the original icosahedron which first rotate 
by 60 in opposite direction with respect to the parallel trian¬ 
gular facets and becomes part of a diamond facet which then 
transforms into a square. 

As first noticed by Macka}f^ an anticuboctahedron 
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could be transformed into an Ih throughout a rotation 
in opposite ways of the two triangular facets perpendic¬ 
ular to the aCo three-fold rotation axis, by 60 degrees 
with reference to each other, along the same axis. The 
transformation is similar with respect to the one of the 
Co, involving six parallelograms, however in this case two 
opposite parallel triangular facets rotate about their nor¬ 
mal and three pairs of abutting triangular facets remain 
unchanged and rotate about the axis that is perpendicu¬ 
lar to their common edge and belongs to the twin plane, 
see Figure 



FIG. 7. Initial, saddle and final configuration in the case 
of AgPt systems. Top row: Ih into Co transformation in 
Ag 92 @Pt 55 . Some surface atoms are removed from the pic¬ 
ture in order to show that the transformation is simultaneous 
in the two chemical species. Bottom row: Dh into Ih trans¬ 
formation in Ag 74 Pt 73 . Ag atoms are in grey and Pt in blue. 
Red lines highlight facets involved in the DSD mechanism. 


The simulation parameters used for the monometallic 
systems are also fruitfully adopted to investigate struc¬ 
tural transformations in AgPt NPs. In the case of a 
Ag 92 @Pt 55 -where stands for a shell/core ordering. 
We would like to note that the results showed in Fig¬ 
ure [prefers to the same system. The pathways followed 
during the structural transitions are the same described 
early. The transition happens simultaneously in the two 
chemical species as shown in the top row of Figure 
Ih configuration is clearly more energetically favourable 
with respect to the Co one. The collective twisting of 
two triangular facets is energetically expensive due to 
the strain increase while the almost barrierless transfor¬ 
mation of Co into Ih geometry is corroborated by Fig¬ 
ure 1^ bottom panel. Roughness of the landscape can 
be bestowed to the sampling of intermediate defected 
structures. Ag 92 @Pt 55 reaches a diffusive regime between 
Ih, Dh and Co configurations before visiting higher en¬ 
ergy conformati ons, a llowing a FES reconstruction ac¬ 
cordingly to Ref b^ l ^Q I. The Ih to Co FES reconstruction 
is shown in Eigurej^ 

A quantitative analysis of the energy barriers defining 
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MPDD 

FIG. 8. Free energy surface reconstruction for an Ag 92 @Pt 55 
run where a diffusive regime in the exploration of Go and 
Ih basins is found. Golor key helps to visualize free energy 
differences, AF, with respect to at a given point in 

the GV space. 


the detected processes was done exploiting our on-the-fly 
analysis tool also when no diffusive regime in the confor¬ 
mation’s space exploration was found. In that scenario, 
the free energy barrier is calculated via the amount of 
gaussians deposited in a simulation starting from spe¬ 
cific initial structures according to the initial point of 
the transition of interest. Results for the monometal¬ 
lic and bimetallic transformations are resumed in Table 
|T] together with the DETPS results in parenthesis. We 
would like to note that the geometric structural transi¬ 
tion pathways, including the saddle configuration, on the 
PES and on the EES are almost identical. 

TABLE I. Based on a TBSMA molecular dynamics, free 
energy barrier by MetaD at 300 K, and the activation en¬ 
ergy calculated by DEBTS, in brackets, for Agi 47 , Pti 47 and 
Ag 92 @Pt 55 , respectively. ”NA” means that that pathway has 
not been observed in MetaD simulations. All the values are 
in eV. 


Mechanism 

Ag 

Pt 

Ag92@Pt55 

Dh ^ Ih 

0.6 (0.45) 

1.9 (2.11) 

1.5 

Ih ^ Dh 

NA (2.71) 

NA (4.55) 

4.3 

Co ^ Ih 

0.4 (0.5) 

1.7 (1.9) 

0.25 

Ih ^ Co 

3.2 (3.2) 

NA (5.08) 

3.6 

aCo —7-Ih 

0.5 (0.64) 

NA (1.87) 

NA 

Ih aCo 

3.0 (3.49) 

NA (5-13) 

NA 


Eree energy barriers heights can be related to the acti¬ 
vation barrier for the diffusion of a Pt and Ag atom over 
Pt and Ag (111) and (100) surfaces. We observe that 
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the energy barrier for processes identifies for the alloyed 
core-shell configurations is in between the one found for 
pure Ag and pure Pt clusters while Pt generally shows 
the highest free energy barriers and Ag the lowest. 

Finally, a 1:1 chemical composition, namely Ag 74 Pt 73 , 
with a core of Pt and a shell Ag-rich, has been consid¬ 
ered. In this case, a Dh transformation into Ih and Ih into 
defected Co for is detected involving both a diamond- 
square-diamond and surface reconstruction mechanisms, 
as depicted in the bottom row of Figure We would 
like to mention that in the case of mixed ordering, fur¬ 
ther studies are needed in order to understand whether 
the combination of these two structural transformation 
mechanisms represents the minimum energy pathway for 
the geometrical phase change of the examined structure. 


IV. CONCLUSION 

We have shown that Metadynamics can be success¬ 
fully adapted to the study of morphological transitions 
in metallic and bimetallic NPs by introducing specific col¬ 
lective variables. These are window function set onto the 
pair distance distribution function, and they refer to the 
stacking fault number and the maximum pair distribu¬ 
tion distance. This set of collective variables has a wide 
range of applicability, from monometallic (e.g. Ag, Pt) 
to bimetallic systems with a small mismatch (e.g. AgPt). 
For nanoalloys with a large mismatch a limitation occurs. 
A broadening of the peaks in the PDDF may result in 
the impossibility of uniquely defining the characteristic 
distances, do in Eq. To recognising transitions occur¬ 
ring between two basins and to estimating the associated 
free energy barrier, even if no systematics diffusive regime 
may be obtained, we have introduced an on-the-fly anal¬ 
ysis based on complementary geometrical and energetic 
order parameters, such as CNA and AEquench- 

We have demonstrated that the solid-solid structural 
transitions among the most common closed shell poly- 
hedra happen following the Lipscomb diamond-square- 
diamond mechanism, which is a collective screw disloca¬ 
tion motion. For the first time, we have shown that the 
Mackay’s description of the inter conversion of an Ih into 
an aCo throughout a DSD mechanism is reproducible. 
We remark that the DSD mechanism appears to be a 
universal pathway as it takes place in monometallic as 
well as in nanoalloys with a small either non negligible 
mismatch, such as in AgArP^ and AgPt. 

Finally, we have demonstrated that the proposed 
MetaD scheme is able to predict free energy barriers in 
a very good agreement with DEPTS activation energy 
barriers. We have found that in bimetallic systems the 
free energy barriers lie roughly in between the values of 
the monometallic cases, although their numerical values 
are not simply rescaled. 
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